Maddie Divine Data Analysis
Madison Pickett
4/25/2022
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ stringr 1.4.0
## ✓ tidyr 1.1.4 ✓ forcats 0.5.1
## ✓ readr 2.1.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(knitr)
library(ggplot2)
library(skimr)
f<- "https://raw.githubusercontent.com/mrpickett26/Helping-Maddie/main/M%20Divine%20VC%20Psychfest%20Data.csv"
d<- read_csv(f, col_names=TRUE)## New names:
## * `Participant ID` -> `Participant ID...1`
## * `Start Date` -> `Start Date...56`
## * `End Date` -> `End Date...57`
## * Progress -> Progress...58
## * `Duration (in seconds)` -> `Duration (in seconds)...59`
## * ...## Rows: 210 Columns: 256## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (26): Semester, Hormones D1 Complete, Hormones D2 Complete, VAS D1 Comp...
## dbl (227): Participant ID...1, Sex, BirthControl, Do you take any chemical c...##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.d <-d %>% dplyr::select(Sex,BirthControl,D2S1_Cortisol,D2S1_Testosterone,D2S1_Progesterone,D2S1_Estradiol,D2S2_Cortisol,D2S2_Testosterone,D2S2_Progesterone,D2S2_Estradiol,D2S3_Cortisol,D2S3_Testosterone,D2S3_Progesterone,D2S3_Estradiol,D2VAS1_Stress,D2VAS1_Shame,D2VAS2_Stress,D2VAS2_Shame,D2VAS3_Stress,D2VAS3_Shame,D2VAS4_Stress,D2VAS4_Shame,D2VAS5_Stress,D2VAS5_Shame,CESD)
d %>% group_by(Sex) %>% summarise_all(~min(.x,na.rm=TRUE))## Warning in min(.x, na.rm = TRUE): no non-missing arguments to min; returning Inf## # A tibble: 2 × 25
## Sex BirthControl D2S1_Cortisol D2S1_Testosterone D2S1_Progesterone
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 Inf 0.103 8.02 0.470
## 2 2 0 0.353 1.29 0.769
## # … with 20 more variables: D2S1_Estradiol <dbl>, D2S2_Cortisol <dbl>,
## # D2S2_Testosterone <dbl>, D2S2_Progesterone <dbl>, D2S2_Estradiol <dbl>,
## # D2S3_Cortisol <dbl>, D2S3_Testosterone <dbl>, D2S3_Progesterone <dbl>,
## # D2S3_Estradiol <dbl>, D2VAS1_Stress <dbl>, D2VAS1_Shame <dbl>,
## # D2VAS2_Stress <dbl>, D2VAS2_Shame <dbl>, D2VAS3_Stress <dbl>,
## # D2VAS3_Shame <dbl>, D2VAS4_Stress <dbl>, D2VAS4_Shame <dbl>,
## # D2VAS5_Stress <dbl>, D2VAS5_Shame <dbl>, CESD <dbl>d<- d %>% mutate(Sex = Sex -1)
d<-d %>% mutate(BirthControl = replace_na(BirthControl,0))
data_v1<- d %>% mutate(
#adjusts minimum for each gender to be 1 by shifting all values over
across(ends_with("Cortisol"),~ifelse(Sex == 1, .+(1-0.1031219), .+(1-0.3446964))),
across(ends_with("Progesterone"),~ifelse(Sex == 1, .+(1- 0.46972), .+(1-0.68095))),
across(ends_with("Estradiol"),~ifelse(Sex == 1, .+(1- 0.31000 ), .+(1-0.31622)))
) %>% mutate(
#takes the log of all the hormones for all 3 periods in D2
across(starts_with("D2S"),~log(.))
)
fn_auc_g = function(s1,s2,s3){
(s1+s2)*40/2 + (s2+s3)*15/2
}
fn_auc_i = function(s1,s2,s3){
(s1+s2)*40/2 + (s2+s3)*15/2 - s1*(55)
}
data_v1 <-data_v1 %>% mutate(D2_Cortisol_AUC_i = fn_auc_i(D2S1_Cortisol,D2S2_Cortisol,D2S3_Cortisol),
D2_Cortisol_AUC_g = fn_auc_g(D2S1_Cortisol,D2S2_Cortisol,D2S3_Cortisol),
D2_Testosterone_AUC_i = fn_auc_i(D2S1_Testosterone,D2S2_Testosterone,D2S3_Testosterone),
D2_Testosterone_AUC_g = fn_auc_g(D2S1_Testosterone,D2S2_Testosterone,D2S3_Testosterone),
D2_Progesterone_AUC_i = fn_auc_i(D2S1_Progesterone,D2S2_Progesterone,D2S3_Progesterone),
D2_Progesterone_AUC_g = fn_auc_g(D2S1_Progesterone,D2S2_Progesterone,D2S3_Progesterone),
D2_Estradiol_AUC_i = fn_auc_i(D2S1_Estradiol,D2S2_Estradiol,D2S3_Estradiol),
D2_Estradiol_AUC_g = fn_auc_g(D2S1_Estradiol,D2S2_Estradiol,D2S3_Estradiol)
)
hist(data_v1$CESD)
hist(data_v1$D2_Cortisol_AUC_i)
hist(data_v1$D2_Cortisol_AUC_g)
hist(data_v1$D2_Progesterone_AUC_i)
hist(data_v1$D2_Progesterone_AUC_g)
hist(data_v1$D2_Testosterone_AUC_i)
hist(data_v1$D2_Estradiol_AUC_i)
hist(data_v1$D2_Estradiol_AUC_g)
# library(GGally)
# ggpairs(data_v1 %>% filter(Sex==1), columns = c("CESD","D2_Cortisol_AUC_i","D2_Testosterone_AUC_i","D2_Estradiol_AUC_i","D2_Progesterone_AUC_i"))
#
# ggpairs(data_v1 %>% filter(Sex==0), columns = c("CESD","D2_Cortisol_AUC_i","D2_Testosterone_AUC_i","D2_Estradiol_AUC_i","D2_Progesterone_AUC_i"))
#
# ggpairs(data_v1 %>% filter(Sex==1), columns = c("CESD","D2_Cortisol_AUC_g","D2_Testosterone_AUC_g","D2_Estradiol_AUC_g","D2_Progesterone_AUC_g"))
# Cesd predicted by sex
m1<- lm(data=data_v1, CESD~Sex)
summary(m1)##
## Call:
## lm(formula = CESD ~ Sex, data = data_v1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.076 -7.076 -2.220 5.924 28.636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.0764 0.7937 44.19 <2e-16 ***
## Sex -1.7128 1.4158 -1.21 0.228
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.525 on 208 degrees of freedom
## Multiple R-squared: 0.006987, Adjusted R-squared: 0.002213
## F-statistic: 1.463 on 1 and 208 DF, p-value: 0.2278data_female<- data_v1%>%filter(Sex==1)
data_male<- data_v1%>% filter(Sex==0)
# p value greater than 0.05 so no significant difference
relCESD<- lm(CESD~BirthControl, data=data_female, na.action = na.exclude)
## Now going to run the full linear model for females, accounting for the interaction effect of birth control for all predictiors, as well as the anova on the linear model
library(broom)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:purrr':
##
## some## The following object is masked from 'package:dplyr':
##
## recodelm_cortisol_bc_i<-lm(data=data_female,CESD~(BirthControl*(D2_Cortisol_AUC_i+D2_Testosterone_AUC_i+D2_Progesterone_AUC_i+D2_Estradiol_AUC_i)))
summary(lm_cortisol_bc_i)##
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Cortisol_AUC_i + D2_Testosterone_AUC_i +
## D2_Progesterone_AUC_i + D2_Estradiol_AUC_i)), data = data_female)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.755 -4.479 -1.179 3.730 22.530
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.78747 1.43076 22.916 < 2e-16 ***
## BirthControl 1.39401 2.69321 0.518 0.60702
## D2_Cortisol_AUC_i 0.06606 0.05693 1.160 0.25136
## D2_Testosterone_AUC_i -0.02858 0.06972 -0.410 0.68360
## D2_Progesterone_AUC_i -0.13369 0.12671 -1.055 0.29646
## D2_Estradiol_AUC_i 0.16577 0.11068 1.498 0.14050
## BirthControl:D2_Cortisol_AUC_i 0.05258 0.10600 0.496 0.62208
## BirthControl:D2_Testosterone_AUC_i 0.10690 0.11422 0.936 0.35382
## BirthControl:D2_Progesterone_AUC_i 0.41069 0.16875 2.434 0.01856 *
## BirthControl:D2_Estradiol_AUC_i -0.66104 0.20976 -3.151 0.00274 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.053 on 50 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.3418, Adjusted R-squared: 0.2234
## F-statistic: 2.885 on 9 and 50 DF, p-value: 0.007966plot(fitted(lm_cortisol_bc_i), residuals(lm_cortisol_bc_i))
hist(residuals(lm_cortisol_bc_i))
qqnorm(residuals(lm_cortisol_bc_i))
qqline(residuals(lm_cortisol_bc_i))
m.aov_i <- Anova(lm_cortisol_bc_i, type = "II")
m.aov_i## Anova Table (Type II tests)
##
## Response: CESD
## Sum Sq Df F value Pr(>F)
## BirthControl 120.5 1 1.8580 0.178968
## D2_Cortisol_AUC_i 185.5 1 2.8610 0.096974 .
## D2_Testosterone_AUC_i 2.7 1 0.0415 0.839424
## D2_Progesterone_AUC_i 88.7 1 1.3672 0.247841
## D2_Estradiol_AUC_i 2.5 1 0.0378 0.846575
## BirthControl:D2_Cortisol_AUC_i 16.0 1 0.2460 0.622079
## BirthControl:D2_Testosterone_AUC_i 56.8 1 0.8759 0.353817
## BirthControl:D2_Progesterone_AUC_i 384.1 1 5.9229 0.018561 *
## BirthControl:D2_Estradiol_AUC_i 644.1 1 9.9317 0.002744 **
## Residuals 3242.7 50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## In this case we see significance in differences in CESD scores from cortisol as a predictor alone (0.9 sig), Progesterone with BC accounting for interaction effects (0.95), and Estradiol accounting for BC interaction effects (0.99)
#Now I will do the same for the AUCg, sorry the variables are fucked
lm_cortisol_bc_g<-lm(data=data_female,CESD~(BirthControl*(D2_Cortisol_AUC_g+D2_Testosterone_AUC_g+D2_Progesterone_AUC_g+D2_Estradiol_AUC_g)))
summary(lm_cortisol_bc_g)##
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Cortisol_AUC_g + D2_Testosterone_AUC_g +
## D2_Progesterone_AUC_g + D2_Estradiol_AUC_g)), data = data_female)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.470 -5.701 -1.551 3.481 23.033
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.25707 10.66338 3.869 0.000317 ***
## BirthControl -33.11539 15.30073 -2.164 0.035243 *
## D2_Cortisol_AUC_g -0.05453 0.04594 -1.187 0.240827
## D2_Testosterone_AUC_g 0.06216 0.05942 1.046 0.300568
## D2_Progesterone_AUC_g -0.02303 0.01732 -1.329 0.189739
## D2_Estradiol_AUC_g -0.07103 0.11128 -0.638 0.526162
## BirthControl:D2_Cortisol_AUC_g 0.07432 0.05491 1.354 0.181964
## BirthControl:D2_Testosterone_AUC_g -0.02272 0.08946 -0.254 0.800563
## BirthControl:D2_Progesterone_AUC_g 0.09045 0.04902 1.845 0.070957 .
## BirthControl:D2_Estradiol_AUC_g 0.33400 0.16912 1.975 0.053818 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.458 on 50 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.274, Adjusted R-squared: 0.1433
## F-statistic: 2.097 on 9 and 50 DF, p-value: 0.04746plot(fitted(lm_cortisol_bc_g), residuals(lm_cortisol_bc_g))
hist(residuals(lm_cortisol_bc_g))
qqnorm(residuals(lm_cortisol_bc_g))
qqline(residuals(lm_cortisol_bc_g))
m.aov_g <- Anova(lm_cortisol_bc_g, type = "II")
m.aov_g## Anova Table (Type II tests)
##
## Response: CESD
## Sum Sq Df F value Pr(>F)
## BirthControl 52.1 1 0.7288 0.39733
## D2_Cortisol_AUC_g 0.7 1 0.0099 0.92097
## D2_Testosterone_AUC_g 98.6 1 1.3776 0.24607
## D2_Progesterone_AUC_g 37.5 1 0.5243 0.47238
## D2_Estradiol_AUC_g 55.1 1 0.7706 0.38422
## BirthControl:D2_Cortisol_AUC_g 131.1 1 1.8321 0.18196
## BirthControl:D2_Testosterone_AUC_g 4.6 1 0.0645 0.80056
## BirthControl:D2_Progesterone_AUC_g 243.5 1 3.4042 0.07096 .
## BirthControl:D2_Estradiol_AUC_g 279.0 1 3.9001 0.05382 .
## Residuals 3576.9 50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##In this case we see significance in differences in Estradiol with BC for CESD score response accounting for interaction effects (0.9), and Progesterone accounting for BC interaction effects (0.99)
##Now moving on to Male AUCg
lm_male_g<-lm(data=data_male,CESD~(D2_Cortisol_AUC_g+D2_Testosterone_AUC_g+D2_Progesterone_AUC_g+D2_Estradiol_AUC_g))
summary(lm_male_g)##
## Call:
## lm(formula = CESD ~ (D2_Cortisol_AUC_g + D2_Testosterone_AUC_g +
## D2_Progesterone_AUC_g + D2_Estradiol_AUC_g), data = data_male)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.679 -6.884 -1.992 5.218 23.004
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.639800 8.466501 4.091 7.95e-05 ***
## D2_Cortisol_AUC_g -0.017297 0.034339 -0.504 0.615
## D2_Testosterone_AUC_g 0.015481 0.036185 0.428 0.670
## D2_Progesterone_AUC_g -0.004648 0.028829 -0.161 0.872
## D2_Estradiol_AUC_g -0.011616 0.071309 -0.163 0.871
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.633 on 116 degrees of freedom
## (23 observations deleted due to missingness)
## Multiple R-squared: 0.003139, Adjusted R-squared: -0.03124
## F-statistic: 0.09131 on 4 and 116 DF, p-value: 0.985plot(fitted(lm_male_g), residuals(lm_male_g))
hist(residuals(lm_male_g))
qqnorm(residuals(lm_male_g))
qqline(residuals(lm_male_g))
male.aov_g <- Anova(lm_male_g, type = "II")
male.aov_g## Anova Table (Type II tests)
##
## Response: CESD
## Sum Sq Df F value Pr(>F)
## D2_Cortisol_AUC_g 23.5 1 0.2537 0.6154
## D2_Testosterone_AUC_g 17.0 1 0.1830 0.6696
## D2_Progesterone_AUC_g 2.4 1 0.0260 0.8722
## D2_Estradiol_AUC_g 2.5 1 0.0265 0.8709
## Residuals 10765.1 116##No reported significant difference in Hormone Levels with respect to CESD scores
#Now moving on to Malue AUCi
lm_male_i<-lm(data=data_male,CESD~(D2_Cortisol_AUC_i+D2_Testosterone_AUC_i+D2_Progesterone_AUC_i+D2_Estradiol_AUC_i))
summary(lm_male_i)##
## Call:
## lm(formula = CESD ~ (D2_Cortisol_AUC_i + D2_Testosterone_AUC_i +
## D2_Progesterone_AUC_i + D2_Estradiol_AUC_i), data = data_male)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.753 -6.793 -2.485 5.806 23.303
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.28293 1.10938 31.804 <2e-16 ***
## D2_Cortisol_AUC_i 0.01094 0.03752 0.292 0.771
## D2_Testosterone_AUC_i 0.03135 0.08092 0.387 0.699
## D2_Progesterone_AUC_i -0.02998 0.05809 -0.516 0.607
## D2_Estradiol_AUC_i -0.06328 0.10349 -0.611 0.542
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.614 on 116 degrees of freedom
## (23 observations deleted due to missingness)
## Multiple R-squared: 0.007213, Adjusted R-squared: -0.02702
## F-statistic: 0.2107 on 4 and 116 DF, p-value: 0.932plot(fitted(lm_male_i), residuals(lm_male_i))
hist(residuals(lm_male_i))
qqnorm(residuals(lm_male_i))
qqline(residuals(lm_male_i))
male.aov_i <- Anova(lm_male_i, type = "II")
male.aov_i## Anova Table (Type II tests)
##
## Response: CESD
## Sum Sq Df F value Pr(>F)
## D2_Cortisol_AUC_i 7.9 1 0.0850 0.7711
## D2_Testosterone_AUC_i 13.9 1 0.1501 0.6992
## D2_Progesterone_AUC_i 24.6 1 0.2664 0.6068
## D2_Estradiol_AUC_i 34.6 1 0.3738 0.5421
## Residuals 10721.1 116##No reported significant difference in Hormone Levels with respect to CESD scores
##Look at stressors with respect to hormones
lm_male_i_stress<-lm(data=data_male,CESD~(D2VAS1_Stress+D2VAS2_Stress+D2VAS3_Stress+D2VAS4_Stress+D2VAS5_Stress))
summary(lm_male_i_stress)##
## Call:
## lm(formula = CESD ~ (D2VAS1_Stress + D2VAS2_Stress + D2VAS3_Stress +
## D2VAS4_Stress + D2VAS5_Stress), data = data_male)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.683 -6.548 -1.816 6.273 21.871
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.53111 1.50543 18.952 <2e-16 ***
## D2VAS1_Stress 0.07294 0.06222 1.172 0.2433
## D2VAS2_Stress 0.08730 0.04349 2.007 0.0468 *
## D2VAS3_Stress 0.03547 0.04704 0.754 0.4522
## D2VAS4_Stress 0.08163 0.11033 0.740 0.4608
## D2VAS5_Stress -0.05930 0.11953 -0.496 0.6207
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.938 on 127 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.1863, Adjusted R-squared: 0.1542
## F-statistic: 5.814 on 5 and 127 DF, p-value: 7.208e-05plot(fitted(lm_male_i_stress), residuals(lm_male_i_stress))
hist(residuals(lm_male_i_stress))
qqnorm(residuals(lm_male_i_stress))
qqline(residuals(lm_male_i_stress))
male.aov_i_stress <- Anova(lm_male_i_stress, type = "II")
male.aov_i_stress## Anova Table (Type II tests)
##
## Response: CESD
## Sum Sq Df F value Pr(>F)
## D2VAS1_Stress 109.8 1 1.3743 0.24327
## D2VAS2_Stress 321.9 1 4.0297 0.04683 *
## D2VAS3_Stress 45.4 1 0.5687 0.45218
## D2VAS4_Stress 43.7 1 0.5474 0.46076
## D2VAS5_Stress 19.7 1 0.2461 0.62069
## Residuals 10145.6 127
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#So for males only timepoint 2 is significant for differences in CESD scores
##Now do this for Shame
lm_male_i_shame<-lm(data=data_male,CESD~(D2VAS1_Shame+D2VAS2_Shame+D2VAS3_Shame+D2VAS4_Shame+D2VAS5_Shame))
summary(lm_male_i_shame)##
## Call:
## lm(formula = CESD ~ (D2VAS1_Shame + D2VAS2_Shame + D2VAS3_Shame +
## D2VAS4_Shame + D2VAS5_Shame), data = data_male)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.888 -7.430 -2.008 6.538 22.252
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.37114 1.03619 31.240 <2e-16 ***
## D2VAS1_Shame 0.06015 0.10131 0.594 0.5538
## D2VAS2_Shame -0.06769 0.06355 -1.065 0.2888
## D2VAS3_Shame 0.11173 0.04996 2.236 0.0271 *
## D2VAS4_Shame 0.09076 0.10034 0.904 0.3674
## D2VAS5_Shame -0.07813 0.10906 -0.716 0.4750
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.366 on 127 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.1064, Adjusted R-squared: 0.07121
## F-statistic: 3.024 on 5 and 127 DF, p-value: 0.01298plot(fitted(lm_male_i_shame), residuals(lm_male_i_shame))
hist(residuals(lm_male_i_shame))
qqnorm(residuals(lm_male_i_shame))
qqline(residuals(lm_male_i_shame))
male.aov_i_shame <- Anova(lm_male_i_shame, type = "II")
male.aov_i_shame## Anova Table (Type II tests)
##
## Response: CESD
## Sum Sq Df F value Pr(>F)
## D2VAS1_Shame 30.9 1 0.3525 0.55376
## D2VAS2_Shame 99.5 1 1.1345 0.28883
## D2VAS3_Shame 438.7 1 5.0007 0.02708 *
## D2VAS4_Shame 71.8 1 0.8181 0.36745
## D2VAS5_Shame 45.0 1 0.5133 0.47505
## Residuals 11141.3 127
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#So for males only timepoint 3 is significant for shame for differences in CESD scores
##Do the same for females
lm_female_i_stress<-lm(data=data_female,CESD~(D2VAS1_Stress+D2VAS2_Stress+D2VAS3_Stress+D2VAS4_Stress+D2VAS5_Stress))
summary(lm_female_i_stress)##
## Call:
## lm(formula = CESD ~ (D2VAS1_Stress + D2VAS2_Stress + D2VAS3_Stress +
## D2VAS4_Stress + D2VAS5_Stress), data = data_female)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.821 -4.527 -1.060 4.895 20.484
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.2051875 2.0130020 14.012 <2e-16 ***
## D2VAS1_Stress -0.0002214 0.0710737 -0.003 0.9975
## D2VAS2_Stress 0.0011820 0.0564798 0.021 0.9834
## D2VAS3_Stress 0.1284010 0.0524936 2.446 0.0181 *
## D2VAS4_Stress -0.0509880 0.0944010 -0.540 0.5916
## D2VAS5_Stress -0.0549777 0.0966732 -0.569 0.5722
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.192 on 49 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.1824, Adjusted R-squared: 0.09896
## F-statistic: 2.186 on 5 and 49 DF, p-value: 0.07081plot(fitted(lm_female_i_stress), residuals(lm_female_i_stress))
hist(residuals(lm_female_i_stress))
qqnorm(residuals(lm_female_i_stress))
qqline(residuals(lm_female_i_stress))
female.aov_i_stress <- Anova(lm_female_i_stress, type = "II")
female.aov_i_stress## Anova Table (Type II tests)
##
## Response: CESD
## Sum Sq Df F value Pr(>F)
## D2VAS1_Stress 0.00 1 0.0000 0.99753
## D2VAS2_Stress 0.02 1 0.0004 0.98339
## D2VAS3_Stress 309.45 1 5.9831 0.01808 *
## D2VAS4_Stress 15.09 1 0.2917 0.59156
## D2VAS5_Stress 16.73 1 0.3234 0.57216
## Residuals 2534.35 49
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#So for females only timepoint 3 (stress) is significant for differences in CESD scores
#shame for females
lm_female_i_shame<-lm(data=data_female,CESD~(D2VAS1_Shame+D2VAS2_Shame+D2VAS3_Shame+D2VAS4_Shame+D2VAS5_Shame))
summary(lm_female_i_shame)##
## Call:
## lm(formula = CESD ~ (D2VAS1_Shame + D2VAS2_Shame + D2VAS3_Shame +
## D2VAS4_Shame + D2VAS5_Shame), data = data_female)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.794 -4.596 -1.396 3.267 22.467
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.6598992 1.3169998 21.762 < 2e-16 ***
## D2VAS1_Shame 0.0735331 0.1513691 0.486 0.629284
## D2VAS2_Shame 0.0002647 0.0742381 0.004 0.997170
## D2VAS3_Shame 0.1479441 0.0406896 3.636 0.000664 ***
## D2VAS4_Shame -0.1725549 0.1129552 -1.528 0.133032
## D2VAS5_Shame 0.0367721 0.1270935 0.289 0.773549
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.869 on 49 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.2541, Adjusted R-squared: 0.1779
## F-statistic: 3.338 on 5 and 49 DF, p-value: 0.01132plot(fitted(lm_female_i_shame), residuals(lm_female_i_shame))
hist(residuals(lm_female_i_shame))
qqnorm(residuals(lm_female_i_shame))
qqline(residuals(lm_female_i_shame))
female.aov_i_shame <- Anova(lm_female_i_shame, type = "II")
female.aov_i_shame## Anova Table (Type II tests)
##
## Response: CESD
## Sum Sq Df F value Pr(>F)
## D2VAS1_Shame 11.14 1 0.2360 0.629284
## D2VAS2_Shame 0.00 1 0.0000 0.997170
## D2VAS3_Shame 623.82 1 13.2199 0.000664 ***
## D2VAS4_Shame 110.12 1 2.3337 0.133032
## D2VAS5_Shame 3.95 1 0.0837 0.773549
## Residuals 2312.21 49
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1###Some things I tried out... didnt like it as much
# Male_Estradiol_g<- data_male$D2_Estradiol_AUC_g
# Male_Test_g<- data_male$D2_Testosterone_AUC_g
# Male_Cort_g<- data_male$D2_Cortisol_AUC_g
# Male_Proges_g<- data_male$D2_Progesterone_AUC_g
# Male_Estradiol_i<- data_male$D2_Estradiol_AUC_i
# Male_Test_i<- data_male$D2_Testosterone_AUC_i
# Male_Cort_i<- data_male$D2_Cortisol_AUC_i
# Male_Proges_i<- data_male$D2_Progesterone_AUC_i
# Male_Estradiol_g_norm <- rnorm(200,mean=mean(Male_Estradiol_g, na.rm=TRUE), sd=sd(Male_Estradiol_g, na.rm=TRUE))
# Male_Test_g_norm <- rnorm(200,mean=mean(Male_Test_g, na.rm=TRUE), sd=sd(Male_Test_g, na.rm=TRUE))
# Male_Cort_g_norm <- rnorm(200,mean=mean(Male_Cort_g, na.rm=TRUE), sd=sd(Male_Cort_g, na.rm=TRUE))
# Male_Proges_g_norm <- rnorm(200,mean=mean(Male_Proges_g, na.rm=TRUE), sd=sd(Male_Proges_g, na.rm=TRUE))
# Male_Estradiol_i_norm <- rnorm(200,mean=mean(Male_Estradiol_i, na.rm=TRUE), sd=sd(Male_Estradiol_g, na.rm=TRUE))
# Male_Test_i_norm <- rnorm(200,mean=mean(Male_Test_i, na.rm=TRUE), sd=sd(Male_Test_g, na.rm=TRUE))
# Male_Cort_i_norm <- rnorm(200,mean=mean(Male_Cort_i, na.rm=TRUE), sd=sd(Male_Cort_g, na.rm=TRUE))
# Male_Proges_i_norm <- rnorm(200,mean=mean(Male_Proges_i, na.rm=TRUE), sd=sd(Male_Proges_g, na.rm=TRUE))
#
# plot(x=data_male$CESD, y=Male_Estradiol_g,Male_Estradiol_i, Male_Test_g,Male_Test_i, Male_Cort_g, Male_Cort_i,Male_Proges_g,Male_Proges_i, na.rm=FALSE)
# main = "Male Hormone Levels",
# at = c(1,2,4,5),
# names = c("AUCg Estradiol", "AUCi Estradiol", "AUCg Testosterone","AUCi Testosterone", "AUCg Cortisol","AUCi Cortisol", "AUCg Progresterone","AUCi Progresterone"),
# las = 2,
# col = c("orange","red"),
# border = "brown",
# horizontal = FALSE,
# notch = TRUE
# data_female<-data_female%>%mutate(
# resid_CESD=residuals(relCESD)
# )
#
# m_null<-lm(data=data_female, resid_CESD~1)
# summary(m_null)
#
# add1(m_null, scope=.~.+D2_Cortisol_AUC_g+D2_Testosterone_AUC_g+D2_Progesterone_AUC_g+D2_Estradiol_AUC_g, test="F")
#
# m2<-update(m_null, formula=.~. +D2_Testosterone_AUC_g)
# add1(m2, scope=.~.+D2_Cortisol_AUC_g+D2_Progesterone_AUC_g+D2_Estradiol_AUC_g, test="F")
#
# m3<-update(m2, formula=.~. +D2_Estradiol_AUC_g)
# add1(m3, scope=.~.+D2_Cortisol_AUC_g+D2_Progesterone_AUC_g, test="F")
# m4<-update(m3, formula=.~. +D2_Progesterone_AUC_g)
# add1(m4, scope=.~.+D2_Cortisol_AUC_g, test="F")
#
#
# lm(data=data_female, resid_CESD~D2_Cortisol_AUC_g+D2_Testosterone_AUC_g+D2_Progesterone_AUC_g+D2_Estradiol_AUC_g)
# summary(m_full)
#
#
# m2<- lm(data=data_female, CESD~BirthControl)
# summary(m2)
#
# m3<- lm(data=data_female, CESD~BirthControl+)
# summary(m3)
#
# m2<- lm(data=data_v1, CESD~ logRange+Migration)
# m3<- lm(data=data_v1, CESD~logRange)
# m4<- lm(data=data_v1, CESD~Migration)
# m5<- lm(data=d, CESD~1)
#
# anova(m1,m2, test="F")R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
#So for females only timepoint 3 (shame) is significant for differences in CESD scores, VERY signifcant (want to emphasize this 0.999)
library(sjPlot)## Registered S3 method overwritten by 'parameters':
## method from
## format.parameters_distribution datawizardlibrary(sjmisc)##
## Attaching package: 'sjmisc'## The following object is masked from 'package:skimr':
##
## to_long## The following object is masked from 'package:purrr':
##
## is_empty## The following object is masked from 'package:tidyr':
##
## replace_na## The following object is masked from 'package:tibble':
##
## add_caselibrary(ggplot2)
theme_set(theme_sjplot())
#Now it only makes sense to look at hormonal female data and discuss it, run the signigicant linear models
#Do this for the female AUCi significant data
lm_cortisol_i_sig<-lm(data=data_female,CESD~(BirthControl*(D2_Cortisol_AUC_i)))
summary(lm_cortisol_i_sig)##
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Cortisol_AUC_i)), data = data_female)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.666 -5.280 -1.236 5.225 24.235
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.78410 1.49281 21.961 <2e-16 ***
## BirthControl 2.03622 2.32926 0.874 0.3857
## D2_Cortisol_AUC_i 0.04766 0.05485 0.869 0.3886
## BirthControl:D2_Cortisol_AUC_i 0.20009 0.09646 2.074 0.0427 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.562 on 56 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.1668, Adjusted R-squared: 0.1222
## F-statistic: 3.737 on 3 and 56 DF, p-value: 0.01611plot(lm_cortisol_i_sig)
plot_model(lm_cortisol_i_sig, type = "diag")## [[1]]
##
## [[2]]## `geom_smooth()` using formula 'y ~ x'
##
## [[3]]
##
## [[4]]## `geom_smooth()` using formula 'y ~ x'
plot_model(lm_cortisol_i_sig, type = "std")
lm_estradiol_i_sig<-lm(data=data_female,CESD~(BirthControl*(D2_Estradiol_AUC_i)))
summary(lm_estradiol_i_sig)##
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Estradiol_AUC_i)), data = data_female)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.001 -6.034 -2.595 4.894 23.816
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.0233 1.4848 22.240 <2e-16 ***
## BirthControl 0.8283 2.5188 0.329 0.7435
## D2_Estradiol_AUC_i 0.1575 0.1195 1.319 0.1926
## BirthControl:D2_Estradiol_AUC_i -0.4383 0.2110 -2.077 0.0424 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.988 on 56 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.0818, Adjusted R-squared: 0.03261
## F-statistic: 1.663 on 3 and 56 DF, p-value: 0.1854plot(lm_estradiol_i_sig)
plot_model(lm_estradiol_i_sig, type = "diag")## [[1]]
##
## [[2]]## `geom_smooth()` using formula 'y ~ x'
##
## [[3]]
##
## [[4]]## `geom_smooth()` using formula 'y ~ x'
plot_model(lm_estradiol_i_sig, type = "std")
lm_proges_i_sig<-lm(data=data_female,CESD~(BirthControl*(D2_Progesterone_AUC_i)))
summary(lm_proges_i_sig)##
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Progesterone_AUC_i)),
## data = data_female)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.506 -6.803 -1.996 5.351 23.929
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.34203 1.46483 22.762 <2e-16 ***
## BirthControl 3.90971 2.49791 1.565 0.1232
## D2_Progesterone_AUC_i -0.07227 0.12008 -0.602 0.5497
## BirthControl:D2_Progesterone_AUC_i 0.31309 0.15362 2.038 0.0463 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.819 on 56 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.116, Adjusted R-squared: 0.06863
## F-statistic: 2.449 on 3 and 56 DF, p-value: 0.07303plot(lm_proges_i_sig)
plot_model(lm_proges_i_sig, type = "diag")## [[1]]
##
## [[2]]## `geom_smooth()` using formula 'y ~ x'
##
## [[3]]
##
## [[4]]## `geom_smooth()` using formula 'y ~ x'
plot_model(lm_proges_i_sig, type = "std")
#Do Do this for the female AUCg significant data
lm_estradiol_g_sig<-lm(data=data_female,CESD~(BirthControl*(D2_Estradiol_AUC_g)))
summary(lm_estradiol_g_sig)##
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Estradiol_AUC_g)), data = data_female)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.413 -6.905 -1.247 4.134 24.811
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.4875 8.9621 4.741 1.5e-05 ***
## BirthControl -27.0616 12.2387 -2.211 0.0311 *
## D2_Estradiol_AUC_g -0.1178 0.1124 -1.048 0.2990
## BirthControl:D2_Estradiol_AUC_g 0.4133 0.1658 2.493 0.0157 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.798 on 56 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.1201, Adjusted R-squared: 0.07297
## F-statistic: 2.548 on 3 and 56 DF, p-value: 0.06497plot(lm_estradiol_g_sig)
plot_model(lm_estradiol_g_sig, type = "diag")## [[1]]
##
## [[2]]## `geom_smooth()` using formula 'y ~ x'
##
## [[3]]
##
## [[4]]## `geom_smooth()` using formula 'y ~ x'
plot_model(lm_estradiol_g_sig, type = "std")
lm_proges_g_sig<-lm(data=data_female,CESD~(BirthControl*(D2_Progesterone_AUC_g)))
summary(lm_proges_g_sig)##
## Call:
## lm(formula = CESD ~ (BirthControl * (D2_Progesterone_AUC_g)),
## data = data_female)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.005 -5.933 -2.128 5.392 23.547
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.10804 3.10247 11.961 < 2e-16 ***
## BirthControl -8.23628 4.18881 -1.966 0.05423 .
## D2_Progesterone_AUC_g -0.02403 0.01706 -1.408 0.16458
## BirthControl:D2_Progesterone_AUC_g 0.12209 0.03809 3.206 0.00223 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.576 on 56 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.164, Adjusted R-squared: 0.1192
## F-statistic: 3.661 on 3 and 56 DF, p-value: 0.0176plot(lm_proges_g_sig)
plot_model(lm_proges_g_sig, type = "diag")## [[1]]
##
## [[2]]## `geom_smooth()` using formula 'y ~ x'
##
## [[3]]
##
## [[4]]## `geom_smooth()` using formula 'y ~ x'
plot_model(lm_proges_g_sig, type = "std")
## Interesting to look at interaction effect of shame and stress of females on CESD at timepoint 3, run a linear model to do that
lm_tress_shame_sig<-lm(data=data_female,CESD~(D2VAS3_Shame*D2VAS3_Stress))
summary(lm_tress_shame_sig)##
## Call:
## lm(formula = CESD ~ (D2VAS3_Shame * D2VAS3_Stress), data = data_female)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.110 -3.855 -1.314 3.252 23.450
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.3822801 1.8244787 16.653 < 2e-16 ***
## D2VAS3_Shame -0.0548903 0.0535437 -1.025 0.31013
## D2VAS3_Stress -0.0450457 0.0496274 -0.908 0.36832
## D2VAS3_Shame:D2VAS3_Stress 0.0027281 0.0009827 2.776 0.00768 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.528 on 51 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.2988, Adjusted R-squared: 0.2576
## F-statistic: 7.245 on 3 and 51 DF, p-value: 0.0003857plot(lm_tress_shame_sig)
plot_model(lm_tress_shame_sig, type = "diag")## [[1]]
##
## [[2]]## `geom_smooth()` using formula 'y ~ x'
##
## [[3]]
##
## [[4]]## `geom_smooth()` using formula 'y ~ x'
plot_model(lm_tress_shame_sig, type = "std")
##Interaction is significant (0.05)Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.